#------------------------------------------------------------------------------
# Import Library
#------------------------------------------------------------------------------

rm( list=ls() )
library( LalRUtils )
libreq(
        tidyverse, data.table, zoo, tictoc, 
        fst, fixest, PanelMatch, patchwork,
        rio, magrittr, janitor, did, panelView, 
        ggiplot, tictoc, binsreg, interflex
       )
set.seed( 42 )
theme_set( lal_plot_theme() )

#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])
#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )

#------------------------------------------------------------------------------






#------------------------------------------------------------------------------
# Generation of graph
#------------------------------------------------------------------------------

vcf = vcf_data[ year>= 1995 ]
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf[, time_fra := year - 2008]
vcf[, D_fra := sch * (year >= 2008)]
vcf[, "t"] <- vcf$year-1995

# %% VCF mining
names(vcf)
vcf_data = copy(vcf)
vcf_data[, `:=`(id_f = as.factor(cellid),
                sty_f = as.factor(styear),
                year_f = as.factor(year),
                D_f    = as.factor(D))]
vcf_data[, distance := min_dist_to_mine * 110]


# %% vcf prep
vcf_estim = vcf_data[cover_1990 >= quantile(cover_1990, 0.5)]
vcf_estim[, mine_dist_q := ntile(min_dist_to_mine, 3)]
setnames(vcf_estim, 'blk', 'block')

vcf_estim[, `:=`(
  D_mine_1 = D * (mine_dist_q == 1),
  D_mine_2 = D * (mine_dist_q == 2),
  D_mine_3 = D * (mine_dist_q == 3))]


# %%
m00_vcf = feols(forest_index ~ D_mine_1 + D_mine_2 + D_mine_3 |
                  cellid + year, cluster = "block", vcf_estim)
# 2wFE + state year FEs
m01_vcf = feols(forest_index ~ D_mine_1 + D_mine_2 + D_mine_3 |
                  cellid + styear, cluster = "block",
                vcf_estim)
m03_vcf = feols(forest_index ~ D_mine_1 + D_mine_2 + D_mine_3 |
                  cellid + cellid[t] + styear, cluster = "block",
                vcf_estim)

vcf_estim[, never_treated := max(D) == 0, cellid]
controls_pre1 = vcf_estim[never_treated == 1 & year < first_pesa_exposure, 
                       mean(forest_index)]
treat_pre1    = vcf_estim[never_treated == 0 & year < first_pesa_exposure, 
                       mean(forest_index)]

# %% GFC
gfc_estim = gfc[pref == 1]
gfc_estim[, mine_dist_q := ntile(min_dist_to_mine, 3)]
gfc_estim[, mine_dist_50_b := ifelse(min_dist_to_mine>5, 1, 0)]
gfc_estim[, `:=`(
  D_mine_1 = D * (mine_dist_q == 1),
  D_mine_2 = D * (mine_dist_q == 2),
  D_mine_3 = D * (mine_dist_q == 3))]


# %%
m00_gfc = feols(def_ha ~ D_mine_1 + D_mine_2 + D_mine_3 |
                  village + year, cluster = "block", gfc_estim)
# 2wFE + state year FEs
m01_gfc = feols(def_ha ~ D_mine_1 + D_mine_2 + D_mine_3 |
                  village + styear, cluster = "block",
                gfc_estim)
m03_gfc = feols(def_ha ~ D_mine_1 + D_mine_2 + D_mine_3 |
                  village + village[t] + styear, cluster = "block",
                gfc_estim)

gfc_estim[, ever_treated := max(D), .(village)]
ctrl_mean_gfc = round(gfc_estim[ever_treated == 0 & 
                               pesa_exposure == 0, 
                             mean(def_ha)], 2)
treat_mean_gfc  = round(gfc_estim[ ever_treated == 1 & 
                               pesa_exposure == 0, 
                               mean(def_ha)], 2)

# %%
mods = list(m00_vcf, m01_vcf, m03_vcf, m00_gfc, m01_gfc, m03_gfc)
# %%
treatmap =c(
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-forest green index",
  "built_index"  = "Non-forest index",
  "cellid"       = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "D_mine_1"     = "Scheduled X PESA X 1st Tercile",
  "D_mine_2"     = "Scheduled X PESA X 2nd Tercile",
  "D_mine_3"     = "Scheduled X PESA X 3rd Tercile",
  "block"        = "Block",
  "blk"          = "Block"
)

# %%
fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")
etable(mods,
       style.tex = style.tex(main = "base", 
                             depvar.title = "", 
                             model.title = "", 
                             var.title = "\\midrule", 
                             slopes.title = "\\midrule \\emph{Time Trends}", 
                             yesNo = c("$\\checkmark$", ""), 
                             signif.code = NA, 
                             tablefoot = FALSE, 
                             line.bottom = "\\midrule \\midrule"),
       signif.code = NA,
       fixef_sizes = T, 
       order = c("cellid", "village", "year" , "styear"),
       fixef_sizes.simplify = F,
       fitstat = c( "n_new", "r2" ), 
       tex = TRUE,
       dict = treatmap, 
       notes = "\\emph{Notes:} Standard errors are clustered at the block level and reported in parentheses.",
       extralines=list(
         "\\midrule \\emph{Summary Statistics}" = c( "", "", "", "", "", ""  ),
         "Mean Pre-Y (Non-Sch)"= c( rep( controls_pre1, 3 ), 
                                    rep( ctrl_mean_gfc, 3 ) ),
         "Mean Pre-Y (Sch )"   = c( rep( treat_pre1, 3 ), 
                                    rep( treat_mean_gfc, 3 ) ),
         "Dataset"     = c( rep( "VCF", 3 ), 
                            rep( "GFC", 3 ) ),
         "Timespan"    = c( rep( "1995-2017", 3 ), 
                            rep( "2001-2017", 3 ) )
       ),
       label = "tab:regs_mining",
       title = glue::glue("Regression estimates decomposed by distance to mines 
                          (ex-ante median cutoff)"),
       file = "appendix_tableA6.tex", replace = TRUE
)

#------------------------------------------------------------------------------
